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DEVELOPMENT AND APPLICATION OF A LOCAL LINEARIZATION 
ALGORITHM FOR THE INTEGRATION OF QUATERNION 
RATE EQUATIONS IN REAL-TIME FLIGHT 
SIMULATION PROBLEMS 

By Lawrence E. Barker, Jr., Roland L. Bowles, 
and Louise H. Williams* 

Langley Research Center 

SUMMARY 

High angular rates encountered in real-time flight simulation problems may require 
a more stable and accurate integration method than the classical methods normally used. 

A study has been made to develop a general local linearization procedure of integrating 
dynamic system equations when using a digital computer in real time. The procedure is 
specifically applied to the integration of the quaternion rate equations. For this applica- 
tion, results are compared to a classical second-order method. The local linearization 
approach is shown to have desirable stability characteristics and gives significant 
improvement in accuracy over the classical second-order integration methods. 

INTRODUCTION 

The quaternion (or Euler parameter) rate equations are widely used for determining 
the orientation of missiles and aircraft in real-time flight simulation problems. These 
equations are commonly integrated by using a second-order integration method developed 
by Adams and Bashforth (AB-2 method). As shown in reference 1, the AB-2 method is 
sufficiently accurate for simulation studies involving aircraft with moderate values of 
angular rates; however, with the advent of high performance aircraft and their subsequent 
simulation, more accurate integration techniques are required. 

The purpose of the present study is to further evaluate integration of the quaternion 
rate equations by the AB-2 method and to develop and evaluate an improved integration 
scheme to use in real-time simulations. Desirable characteristics for this improved 
algorithm are: (1) stability and accuracy over a large range of angular rates, (2) one- 

pass algorithm which could be used at large step sizes, and (3) computer timing and 
memory requirements comparable to the classical methods commonly used. 

*Electronic Associates, Inc. 



A general algorithm based on a local linearization procedure is developed and is 
then applied to the digital integration in real time of the quaternion rate equations. The 
resulting algorithm is referred to as the LL algorithm. The description of the LL imple- 
mentation on a digital computer is presented in the form of a program flow chart and 
FORTRAN source listing. Also included is a simplified version of the LL algorithm for 
those users with limited computer resources. 

The LL algorithm is compared with the AB-2 method both analytically and experi- 
mentally. This comparison includes stability and error analysis. Included for the exper- 
imental study are both analytical inputs and inputs taped from an actual piloted run. In 
addition, both methods (LL and AB-2) are compared with a variable step seventh-order 
Runge-Kutta method used as a high-quality approximation to the exact solution. 

SYMBOLS 


A\B\C\C^,C^ 
C it D',H,G,J,K 


coefficients in LL quaternion algorithm, where i = 1, . . .,4 


A(t) time-varying matrix; for quaternions, a skew-symmetric matrix which relates 

components of a vector to their derivatives due to referencing a rotating 
reference frame, radians/second 

Ap amplification constant for roll rate p, radians/second 

a i quaternions (Euler parameters), where i = 1, . . ? 4 

B control matrix of dynamic system 

C transformation matrix relating vehicle body axes to inertial axes 

Cjj matrix elements of C, where i,j = 1, 2, 3 

A kh 

e discrete form of transition matrix 


F(X,t) n-dimensional vector composed of general nonlinear time-varying functions 
of state vector X and time t, per second 


9F 

ax 


Jacobian matrix of vector F with respect to vector 


X, 


8F 8F 

aXi’ ax^’ ' ’ 


8F 

9X n 


2 



h 

I 

i 


*i 

k,n 

M 

N 

0(h n ) = (p 
p k ? Qk 


P,q,r 


t 

*k 

u(t) 

X(t) 

X N 

6t 

6X(t) 

X 


integration interval size, seconds 
identity matrix 
imaginary unit, ^1 

constant complex vectors, where i = 1, 2 
constant real vectors, where i = 1, 2 
integer constants 

matrix which relates quaternions at t - ^+1 to those at t = t^ 
norm of X 

where |<f>| = k| h n | and K > 0 

coefficients in general LL algorithm for solution of nonlinear time-varying 
dynamical system 

components of angular velocity vector uJ, radians/second 
time, seconds 

value of t at t = kh, where k = 0, 1, 2, . . seconds 

control vector (forcing function of time-varying system) 

n-dimensional vector representing system states of a dynamic system 

X whose elements have been normalized 

small perturbation in t about t^ 

small perturbation in X about Xjj 

complex number, X = R(X) + I(X)i, radians/second 
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p 


k 


*P,0,<P 


U) 


( j ) 


c 


U)' 


D 


a; 


Where R 

R k 

R 

R AB2 


r ll 

AR 


e R 

^max 

R o ,R(0) 


Where R 


R 

R 

R 


k 

-1 

T 


parameter in LL quaternion algorithm, ('o>j c hj/2, radians 

Euler angles relating vehicle body axes to inertial axes system, radians 

magnitude of u>, radians/second 

actual computed frequency of quaternions, radians/second 

desired frequency of quaternions, cn/2, radians/second 

angular velocity vector relative to inertial axes system, radians /second 

is any arbitrary scalar or vector variable: 

value of R at t = t^ 

derivative of R with respect to time 

R evaluated by AB-2 method 

R evaluated by LL algorithm 

deviation in R from some reference 

error in R from some reference 

maximum value of R 

initial value of R 

is any arbitrary matrix: 

value of R at t = t^ 

inverse of R 

transpose of R 
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ABBREVIATIONS 


AB-2 Adams -Bashforth second-order integration 

LL local linearization 

RK-7 Runge-Kutta seventh-order integration 

ANALYSIS 

Statement of Problem 

An aircraft may be considered a moving coordinate system. From the angular 
velocity components (p,q,r), an integration of angular rate to determine angular position 
is done. The quaternion rate equations are normally used in describing this orientation 
(refs. 2 and 3). In some problems such as the simulation of two aircraft, the relative 
geometry computations and the integration of the equations of motion, as well as the above 
mentioned angular rates, require a large number of arithmetic operations per integrating 
step. These operations are so numerous that an integration method is needed that 
requires only one evaluation of the derivatives per integration step. The term "time 
critical" will be used to define problems with this characteristic. A method which does 
n derivative evaluations per step is defined as an n-pass integration scheme. Higher 
order multistep methods could be used to increase the accuracy and maintain one evalua- 
tion per step, but these methods have other difficulties such as determining the required 
starting values. At the NASA Langley Research Center an interval size of h = 1/32 sec 
is normally used for real-time flight simulation problems. The integration is typically 
done with a one-pass integration method such as second-order AB-2. For high angular 
rates encountered in high-performance aircraft, particularly as subsystems of an air-to- 
air combat simulation, the computation is not sufficiently accurate using this method. 

The purpose of the present study is to develop an improved algorithm for the integration 
of the quaternion rate equations on a digital computer in real time. For this report, 
angular rates up to 10 rad/sec are considered. 

Figure 1 shows how the quaternion equations fit into a simulation problem. Note in 
this figure the direction cosines c^j are computed algebraically from the quaternions. 

The angular components p, q, and r of the spin vector to, available from the equations 
of motion, are taken as inputs to the quaternion block which is treated as the system in this 
study. Also shown are the initialization and attitude angle readout blocks. The remaining 
block in the figure is the normalization block. This block is usually added to normalize 
the quaternions after integration of the quaternion rate equations (appendix A). Any 
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closed-loop effects through the equations of motion or any visual or other feedback effects 
due to a pilot are not considered in this report. 

The quaternion equations can be written as the following matrix differential 
equation: 


X(t) = A(t) X(t) 


( 1 ) 


where 

X(0) = x 0 


0 



-r -q -p 

0 -p q 

p 0 -r 


p -q r 0 



( 2 ) 


a l 

_ a 2 

x = (a 

a 3 

a 4 

_ J 

(Note that aj, ag, a^, and a^ of this report are equivalent to ag, ag, ag, and a^, 
respectively, of ref. 3.) From the eigenvalues of the above system, it is seen that the 
natural frequency of this system is co/2, where 

< j ? - p^ + q^ + r 2 


Other methods of calculating attitudes of a rigid body such as direction cosines (refs. 2 
and 3) result in a frequency of w. The reduced frequency is an additional advantage in 
using the quaternion equations. 

For the matrix A(t) defined in equation (2), matrix equation (1) is called the Euler 
parameter equation where a p . . ., a 4 from equation (3) are called Euler parameters 
or quaternions. The normalization condition for the quaternions is defined as 
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(4) 


N 


2 



= 1 


where N is the norm. The term "norm in" will be used to define the system described 
by equation (1) and satisfying equation (4). When the system is not required to satisfy 
equation (4), it will be referred to as "norm out." The effect of normalization upon 
accuracy of computation is discussed later. This report is concerned with the integra- 
tion of equation (1) in a one-pass system where A(t) defines the quaternion dynamics. 
The investigation is performed for the classical second-order Adams-Bashforth predictor 
integration method (AB-2 method) and for the new algorithm (LL algorithm) which was 
developed to maintain accuracy at high angular rates and at increased step sizes. 


AB-2 Method 

Before developing the LL algorithm, a look at the solution of equation (1), using the 
AB-2 method with Euler integration as the starter formula, is in order (ref. 1). The vec- 
tor AB-2 equation is written 


X k+1 " X k + |( 3X k ‘ X k-l) 


Applying equation (5) to equation (1) gives 


X, + -(3 A, X, - A. .X,. 


k+r A k + 2 J Vk - "k-i^k-i 


(5) 


( 6 ) 


where 


A k = A(t k ) 


Where w is a constant vector A = A k , by taking the z-transform of equation (6) 
and solving the resulting characteristic equation as in reference 1, the approximate solu- 
tion assuming ho> « 1 can be written as 


X(t) 


4h 3 a)4 

256 


.to 5,23). ./<*> 5 .2 3) , 

i — +— h or t -i — +— h t 

- \2 96 I 7 o 2 96 

K^e x + I^e ' 


(7) 


where and are constant complex vectors that depend on initial conditions. The 

application of Euler's formula yields the real solution 
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X(t) = 


4h3c£^ t 
, 256 


k; 


cos 


(f 


_5_ 

96 


h 2 u> 3 


t + 




sin/— + A h 2 ^ 3 ^ 
\2 96 ) 


( 8 ) 


where Kj and are constant real vectors. Using initial conditions for equation (1), 
equation (8) becomes 


X(t) = e 


4h 3 uA 

256 


I cos 


\2 96 


h^ct)3u + 


4h^o?^ 

256 


— + — h 
2 96 


2 a,3 


sin| ,M h2 


w-» t 


X(0) 


which can be approximated by 


X(t) 


4h 3 o) 4 t 
„ 256 1 


I cos/— + — h 3 cu 3 )t 
\2 96 



(9) 


The solution for the norm squared is found from equation (9) to be approximately (ref. 1) 


N 


2 



( 10 ) 


This equation will be used later in discussing the error in the norm. 


General Local Linearization Algorithm 

A general algorithm for the solution of a nonlinear system is now developed by the 
method of local linearization. This approach is taken to produce a single-step one-pass 
integration algorithm with the intent of applying it to the quaternion rate equations where 
classical schemes fail for large rates. As mentioned previously, a one-pass scheme is 
very important for time-critical problems such as large digital simulations in real time. 
The algorithm is derived by applying the method of perturbation to the differential equa- 
tions (not states) and then solving exactly the resulting differential equations after dropping 
higher order terms in the perturbation. 

Let the nth order ordinary differential equation describing a general nonlinear time- 
varying dynamical system be expressed in state variable form as 


8 



X(t) - F(X,t) 


(ID 


where X(t) is the n-dimensional vector representing the system states and F is the 
n-dimensional vector composed of general nonlinear time-varying functions of the state 
vector X(t) and time t. Let 6X(t) be a small perturbation about X k defined as 

£>X(t) = X(t) - X k 


where t k = t = t k and let 


6t be a small perturbation in time about t k defined as 


6t = t - t k 


where t k = t = t k+ ^- Therefore 


6X = X = F(x k ,t k ) 


H(S k> t k )6X(t) + ^ 

0X V j 0t 


: (%tk) 


6t 


(12) 


where 

(Vk) 


F(X,t) has been expressed as a Taylor series in 
When equation (12) is integrated, it becomes 


X(t) 


and t about the point 



where 


A k -A 





h = *k + l " 


or 


x k + i - x k - p k F ( x k- i k) + «k f( x k.‘k) - 4 3 ) 


(13) 
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where 


P 


k 





- I - hA,, 


Equation (13) is a general algorithm for the solution of the nonlinear problem of equa- 
tion (11) and is related to the exponential method of reference 4. Note that this is a 
single-step solution to the nonlinear problem with local truncation error being o(h 3 ). 

Now, assuming 

F(X,t) = A(t) X(t) + Bu(t) 

(a linear time-varying system with forcing function u), equation (13) becomes 

*k+i = e k X k + Q k A A + P k Ba k + < 14 ) 

For B = 0 (a time -varying system with no forcing function), equation (14) becomes 


X k + 1 ' * «k\ X k »5) 

Equation (15) is the particular case that is applied in this report to derive an integration 
algorithm for the quaternion rate equations. 


LL Algorithm 

Application of equation (15) to the quaternion differential equations results in the 
LL algorithm. For the matrix series definition of e Ak , it can be shown that 



= I cos 


2 


u) k h 

+ — - sin — 
w k 2 


(16) 


since 


- 



I 


( 17 ) 
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where 




2 - 2 + r 2 


+ \ p k +q k +1 k 

Substituting equation (16) into equation (15) and using equation (17) results in 


X k+1 " C l X k + C 2 A k X k + C 3 A k X k + C 4 A k A k X k 


( 18 ) 


where the coefficients are defined by 


C x = cos p k 
2 sin 


Cq = 

l (j) 


C 3 = o\ ‘ C ° S P k 
W k 


and 


C A = — — (h — — sin p, 

4 u,. 2 V “k k 


“k h 


Pk 2 


With matrix A defined by equation (2), equation (18) can be written 


X k + 1 = M k X k 


(19) 


where 


M,, 


H 

-G 

J 

K 


G 

H 

K 

-J 


-J 

-K 

H 

-G 


-K| 

J 

G 

H 
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and where the scalar elements of M, are 

k 


H = Cj + C 4 A’ 

(20) 

G = " C 2 r k " C 3 r k + C 4 B 

(21) 

J^k + ^k-V’ 

(22) 

K = C 2 p k + C 3 p k " C 4 D 

(23) 


A ._ PkPk + %% + _ Vk 

4 

B -.^kVVk 

4 

c . _ Vk - Pk^k 
4 

D . _ q k*k - qk r k 
4 


C 3 " ' cos p k) 

k 

Equation (19) constitutes the LL algorithm. 

The FORTRAN subroutine (QUAT) incorporating the LL algorithm plus normaliza- 
tion is given in appendix B. As mentioned before, this is a single-step algorithm. (See 
fig. 2 to compare its implementation with that of the AB-2 method.) Note that in addition 
to using the spin vector w, this algorithm uses the angular acceleration Z 5, which is 
always available from the equations of motion. This adds to its accuracy. If stringent 
requirements are necessary, a simplified version of this algorithm requiring less time 
and memory may be used. See appendix C. 
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Assuming oj is a constant vector, equation (15) becomes 


X - e^X 
A k+1 " e A k 


which is the discrete form of 
X(t) = e At X(0) 


( 24 ) 


or 


X(0) (25) 

where equation (16) was substituted into equation (24). In contrast to the AB-2 method, 
equation (9), equation (25) gives the exact solution for the LL algorithm. Refer to appen- 
dix D for further stability characteristics of the LL algorithm. 

Computational Considerations 

The LL algorithm developed in this report, along with the AB-2 method, has been 
implemented on a CDC 6600 computer at Langley Research Center and has been pro- 
gramed to operate in a real-time mode. Most runs were made at a step size of 
h = 1/32 sec, which is the normal step size that is used for real-time problems. The 
implementation of the LL algorithm is simple. The flow chart of figure 3 depicts the flow 
of the real-time program and shows the placement of the LL quaternion subroutine (sub- 
routine QUAT). It should also be noted that the body rate input variables p, q, and r, 
and their derivatives p, q, and r, will be stored by subroutine SAVE during the first 
pass of a multipass integration scheme. For a single-pass scheme this storage may be 
omitted. 

Additional compute time for calculating the exponential matrix function of equa- 
tion (15) by summing a series or by other numerical methods did not exist for the partic- 
ular application of this report. The skew-symmetric nature of the quaternion system of 
equation (1) enabled a closed analytical expression to be written. However, the LL algo- 
rithm does have slightly larger core and timing requirements than the AB-2 method. The 
additional requirements were approximately 88 locations and less than 1 percent of the 
compute time per step size (h = 1/32 sec). These requirements are negligible for the 
typical real-time flight simulation problems; these requirements are even more insignifi- 
cant when the accuracy gained and the possibility of significantly increasing the step size 
are considered. 


X(t) = 


I cos(- t) + — sin/— t 

cu 
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RESULTS AND DISCUSSION 


Experiment 

Two methods of computing the quaternions (AB-2 and LL) were incorporated into a 
real-time program. These methods were evaluated simultaneously and the results were 
compared with each other and with the results obtained independently from a high quality 
seventh-order variable-step Runge-Kutta integration routine. Most runs were made at 
h - 1/32 sec. As mentioned previously, the quaternion equations are the system or plant 
of concern for this experiment and require as inputs the angular velocity components p, 
q, and r of the spin vector tu. (The LL algorithm, in addition, requires the angular 
acceleration components p, q, and r.) The results will be given for four different input 
types, progressing from the simpler constant inputs to an actual taped piloted run. These 
inputs are: (1) constant rates, (2) sinusoidal rates, (3) pulse sinusoidal rates, and 

(4) piloted inputs. 

Figure 1 shows the variables of interest in this study. In addition to the errors in 
the quaternions, the errors encountered in typical direction cosines and in the Euler 
angles ip, 6, and 6 are given. The error in <b is of special interest due to the vehi- 
cle roll rate p being the main contributor to the angular velocity vector cu. These 
errors would be important in the overall flight simulation problem since they are sent to 
the simulator and fed back to the equations of motion through pilot control inputs. Though 
important, the effects of the errors on the complete simulation problem will not be given 
in this report. The effect of "norm out” will be given for the constant rate input for both 
the LL and AB-2 methods for comparison (see ref. 1). Results for the other inputs will 
be for "norm in” unless stated otherwise. 

Constant inputs .- The simpler constant rate inputs were chosen since they permitted 
an analytical solution (eq. (9)) and allow comparison with the results of reference 1. In 
addition, this case allows the LL algorithm to be used as a reference since it is exact, 
neglecting computer round-off errors. 

Figure 4(a) shows for the AB-2 method the percent error e N 2 , from its desired 
value of 1.0. (See eq. (4).) This error is plotted as a function of time for the constant 
values p = q = r = l/^3 rad/sec and h = 1/16 sec. Normally the Euler integrating 
starting formula is used to start the AB-2 method. The starting error has been sub- 
tracted as in reference 1. Figure 4(b) shows the norm squared error for h = 1/32 sec 
for p = 1, 3, and 5 rad/sec (q = r = 0). As shown in the figure for p = 1, there is only 
negligible error after 120 sec, but as p is increased, larger errors are encountered at 
shorter times. Thus the curves indicate the AB-2 method may no longer be applicable 
for these larger roll rates. The predicted errors shown were calculated from equa- 
tion (10). The errors in the norm can be eliminated by the addition of normalization to 
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the resulting quaternions. As stated previously, the LL algorithm has no error in norm 
for a constant angular velocity. It might be concluded that the LL algorithm does not 
require normalization. However, for the general case, cu varying, this normalization 
is shown to add a few digits of accuracy. 

See figures 5(a) and 5(b) for two time-history segments of the quaternion a^ for a 
constant roll-rate input, p = 10 rad/sec (q = r = 0). The upper time-history plot in 
each figure represents a^ resulting from the LL algorithm, and the lower represents 
a^ resulting from the AB-2 method. Note, in figure 5(a) the divergence of aj for AB-2 
which could have been predicted from the results shown in figure 4(b). When normaliza- 
tion is added, figure 5(b) shows the error in norm has been eliminated for AB-2. 

Figure 6 shows the error in computed frequency of the quaternions e for AB-2 

plotted against the correct or desired frequency This frequency error was invariant 

of normalization. The solid curve shows the error as predicted by equation (9). The 
errors for LL were essentially zero and are thus not plotted. 


Figure 7 shows for AB-2 the error in roll angle versus time for constant values of 
roll rate p and h = 1/32 sec. Even for values of p where is small, the direc- 
tion cosines could still have an effect in the inertial velocities or relative geometry. (See 
fig. 1.) 


Figure 8 shows how the errors in the direction cosines vary with time. Since the 
errors for LL are essentially zero for constant rate inputs, LL is used as a reference in 
calculating these errors for AB-2. The beat oscillation with period of 62.8 sec is due to 
the phase shift term given in equation (9). Not all elements C- (where i,j = 1, 2, 3) 
of the direction cosine matrix C vary since q = r = 0 for this case. The last two time 
histories show typical direction cosines, C 22 and C 32‘ 

Sinusoidal inputs .- The following sinusoidal inputs were chosen to provide a drastic 
time-varying case, drastic in the sense that w is approximately 10 rad/sec: 


p = 10 sin 0.5t 


q = r = 2 sin t 

Figure 9(a) shows typical time histories of the quaternions a^, . . ., a^, while fig- 
ure 9(b) gives the Euler angles, \p, 6, and </>. A seventh-order Runge-Kutta solution 

(RK-7) was used as an independent check. The continuous curves for the AB-2 method 
and the LL algorithm are denoted by subscripts AB-2 and LL, respectively. Figure 9(c) 
shows the differences in solution between AB-2 and LL for the direction cosines. Also 
shown are two typical direction cosines Cg 2 and Cgg- Figure 9(d) gives the difference 
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in solution for the Euler angles. Also shown are the angular rate inputs p, q, and r. 

The Euler angle errors are given in percent of one revolution (360°). 

Figures 9(a) and 9(b) are plotted on a scale such that any significant error may not 
be readily apparent. However, if reference is made to table I, which gives typical errors 
for t = 58, 59, and 60 sec, significant errors will be noted. (For example, at t = 58 sec 
the AB-2 method shows an error of approximately 14.5°.) The errors are shown for 
AB-2 with "norm in" and for LL with both "norm in" and "norm out." This table clearly 
points out the superiority of the LL algorithm over the AB-2 method. 

The extent of the nonorthogonality of the direction cosine matrix resulting from the 
quaternion integration errors is shown for this time-varying case. The criterion used to 
evaluate the orthogonality of the direction cosines is to multiply the computed matrix C 
of direction cosines by its transpose (appendix A). For the above inputs, the ele- 

ments of the matrix CC T is compared with the corresponding elements of a unity matrix. 
At the end of a 60-second run for AB-2 and "norm out," the product of the computed matrix 
of direction cosines with its transpose produced a matrix with numbers along the main 
diagonal that differed from unity by about 8 x 10“*; for the same run for LL and "norm 
out" the difference was 2 x 10“^. The off-diagonal numbers which should have been zero 
contained only the round-off error of the computer for both methods (4 x 10-15). At the 
end of a 60-second run for "norm in,” the matrix CC^, neglecting round-off error, was 
equivalent to the unity matrix for both methods, as would be expected. 

Sinusoidal pulse inputs .- Results for two sinusoidal pulse cases were chosen to 
provide time-varying inputs with characteristics similiar to those of an actual piloted 
roll maneuver. The roll rate was chosen to give roll in one direction only to show the 
effect of cumulative error in roll angle. 

For the first case, a coning run, where, for positive values only, 
p = 5 sin 0.25t 


with 


q = 0.25 cos 12t 
r = 0.25 sin 12t 

was considered. Figures 10(a) and 10(b) show the quaternions and Euler angles after 
t = 60 sec. Figures 10(c) and 10(d) show how the errors increase with time for both the 
direction cosines and Euler angles. The cumulative error in <f> can easily be seen. 
Typical errors are given in table II. 
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In aircraft studies, since the vehicle roll rate p is the main contributer to the 
angular velocity vector To, a design chart was developed for the second case, a pure roll 
maneuver using AB-2, where, for positive values only, 

p = Ap sin cut 


with 


q = r = 0 


and 


Pmax = A p = 10 rad/sec 
P max = Apcu S 10 rad/sec 2 

Figure 11 is a design chart showing curves of constant error in roll angle <p for two 
consecutive sinusoidal pulses in roll rate p using AB-2. This error is percent full 
scale where full scale is one revolution (4> = 360°). For example, given for a particular 
problem that p> = 3.5 rad/sec 2 and p = 7 rad/sec, the error in roll angle is 
approximately 3 percent or 10.8°. From this chart, one can see that the main contribu- 
tion to error is large values of the angular rate. Even for small accelerations, very 
large errors will result for large angular rates. The errors for this chart were obtained 
using LL as the reference. This was done for programing simplicity since the maximum 
error for LL for the ranges of p and p shown in the chart is 0.1 percent (or 0.36°) for 
Pmax = ra d/sec; thus LL is an adequate reference. Similiar design charts were devel- 
oped for the second-order Adams-Moulton and Runge-Kutta integration routines, which, 
even with their additional pass, did not eliminate the errors at the higher roll rates. This 
could have been predicted from the stability chart in figure 12. 

Piloted run .- The final case used as inputs the angular rates and accelerations from 
an actual piloted run. This was not a very violent run, having a peak roll rate of only 
5 rad/sec for a short duration of time. Figures 13(a) and 13(b), depicting errors in direc- 
tion cosines and Euler angles for AB-2 for a time segment of 130 sec, show p to be the 
dominant input, the criterion used in selecting the previous inputs. Typical errors for 
both AB-2 and LL are given in table III. 
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Effects of Increasing Step Size 


Of great importance in real-time digital flight simulation would be the increase in 
step size without significant loss in accuracy. For the drastic sinusoidal case, fig- 
ures 14(a) and 14(b) show the results for h = 1/16 sec. Similar results for h = 1/32 sec 
were shown in figures 9(a) and 9(b). Typical errors are given in tables I, II, and III for 
the three time-varying cases covered in this report. On a design chart (similiar to 
fig- 11) for the LL algorithm for h = 1/16 sec, the maximum error in <f> was 0.4 per- 
cent (1.4°). These results show the LL algorithm may be used with a large step size to 
significantly decrease the total computing time for a solution and still give results which 
are an order of magnitude more accurate than the AB-2 method using half the step size. 
For systems other than the quaternions, the F(X,t) of equation (1) may be a more com- 
plex vector function and the matrix exponential function e^^ may not be computed 
directly as in this report; however , the possibility of using a larger step size could still 
outweigh the increase in computation time necessary to compute the exponential function. 

CONCLUDING REMARKS 


A general algorithm for the solution of nonlinear systems is developed by the method 
of local linearization. Applying this general formula, a one-pass algorithm has been 
developed for the integration of the quaternion differential equations used in real-time 
digital simulations for six-degree-of-freedom problems. This algorithm has local trun- 
cation error of order of the cube of the step size, making it superior to those classical 
second-order integration schemes which have error of order of the square of the step 
size. Its superiority, in terms of both accuracy and stability, has been demonstrated for 
large angular rates. In fact, it can be used with larger step sizes and still retain accuracy 
and stability. Being a one-pass algorithm, its application to time-critical flight simulation 
problems is possible and sometimes necessary. 

A FORTRAN subroutine has been written to implement this algorithm and its mech- 
anization is very simple. It had a slight increase in timing and memory requirements 
over the classical second-order integration method used for comparison in this report; 
however, this increase was insignificant considering the accuracy gained. If stringent 
requirements are necessary, a simplified version of this algorithm requiring less time 
and memory may be used. Where large angular rates are anticipated or runs of long 
duration are to be made, the algorithm should be mechanized. Even without these condi- 
tions, the algorithm is as efficient and much more accurate than the classical methods 
previously used. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., September 26, 1973. 
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APPENDIX A 


NORMALIZATION 


Numerical integration of the direction cosine rate equations produces errors in the 
resulting direction cosines which in turn cause the resolved components of a given vector 
to be nonorthogonal. Normally, constraint equations on the direction cosines are added 
to improve the orthogonality of the transformation between body and inertial coordinates 
(refs. 5 to 7). In the present report the quaternion rate equations are integrated and the 
direction cosines are then computed algebraically from the resulting quaternion states 
a l> • - a^. This integration may produce errors in the resultant quaternions which 

cause the components of the quaternion vector X of equation (3) to violate the normality 
equation (4). The ultimate result of this violation is the nonorthogonality of the transfor- 
mation matrix mentioned above. When quaternion rate equations are used, steepest 
descent or other methods are employed to satisfy the normality constraint of equation (4) 


(refs. 2 and 3). In the present report the quaternions aj, . . . , a^ are normalized after 
integration of the rate equations. The implementation is accomplished by dividing the 

/ T \ 1 / 2 

quaternions by the norm (X XI 




H 

N 



This mechanization is shown in figure 1. The effectiveness of the constraint can be found 
by multiplying the computed direction cosine matrix by its transpose. For an orthogonal 
transformation, a unit matrix is the correct result of this multiplication. 
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APPENDIX B 


LL QUATERNION SUBROUTINES 


c 

C * * * * 

c * * * * 


OF T H c QijATrRNIOM f^RAMp T P RS 


SUBROUT I NE OUAT ( I ) 

SUBROUT I NE FOR f'.'vALUAT 1 ON 
USING LOCAL LINEARIZAT I ON 
COMMON 

X /I NT COMM/ T , H 

X I SCHEME * DtwlNT(?i 

X ZBKO P / ABM?) ♦ AR2(2) 

X /TOUAT / PK(?) , QK ( 2 ) 

X G00T;<<2) • RD0TK(2> 

DIMENSION OMr.GA;<(2)i OMEGA K 2 ( 2 ), RHOK ( 2 ) « 
X C4 (2 ) • AK (?) ♦ => < ( ? ) , 

X HK (2 ) * SRHOK { 2 ) » CPHOK (? ) 

RFAl JK(2)» KK(2) 

EQUIVALENTS ( r 1 « C RHOK ) 

DATA OM2ZERO / l.E-6 / 

OMFGAK2 ( I ) = P< ( I ) -*P< { I ) + QK ( I ) #Q< ( I ) 4 

IF ( P m E G A K ? ( I ) .{_ T • OM2ZERO ) OMPGAK? ( I ) - 


I NT 
) 

A fp ( 2 ) 

RE ( ? ) 

■1 ) « Cl(2) 

* r< ( p ) 


► NE Q 

» AR4 ( 2 ) 

1 PDOTK ( 2 ) • 

C2P (2 ) iC 3P ( P ) < 
DK ( 2 ) , 0 ^( 2 ) , 


(I > *RK ( I ) 

M2 ZERO 


OMr G A< ( I ) 
RHOK ( I ) 
SRHOK ( I ) 
CRHOK ( I ) 

c;?p ( 1 ) 

C3 P ( I ) 


SORT ( OMEGA <2 ( I) ) 

H#OMEGA< (I ) * # r 
3 I N (RHOK ( 1 ) ) 

COS (RHOK ( I ) ) 

SRHOK ( I ) /OMFGAK ( I ) 

?•*<!• - CRHOK < I ) )/0M^GAK2 ( I ) 


C4 

( I ) 

= 

4.4 (H - 2 • *C?P ( I )) /OMFGAK? ( 

1 ) 


^ < 

( I ) 

= 

- .23* 

(PK ( I ) *PQQTK ( I ) 

4 OK ( I ) 

*Qnn TK ( I ) 

4 RK ( [ ) *RDOTK ( 

r K 

( I ) 

= 

• 2 3 * 

( P< ( I ) *GROTK ( I ) 

- OK ( I ) 

*pro r< t 1 > 

) 

CK 

( T ) 

= 

• 23* 

( RK ( I > *POO TK ( I ) 

- PK ( M 

*RH 0 t K < I ) 

) 

DK 

( I ) 

= 

* 2 3 * 

{ Q;< < I ) *R'r>OTK ( I ) 

- R< ( I ) 

*QQO TK < I ) 

) 

HK 

( I ) 

= 

Cl ( I 

) 4 C4 < I > * AK ( I > 




OK 

( I ) 

- 

-C2R ( 

I ) *RK ( I ) - C3P ( 1 

1 ) *ROOTK 

( I ) 4 C 4 ( 

I ) *BK { I ) 

JK 

( r ) 

= 

C2P ( 

I )*OK ( I ) 4 r 3P ( ] 

t ) *QDOTK ( I ) - 04 ( 

I ) *CK ( I ) 

KK 

( 1 ) 

= 

C2P ( 

I ) *PK ( I ) 4 C 3P ( i 

I ) *PDOTK 

( 1 ) - C 4 ( 

I ) *DK ( I > 


ATI 
A T 2 
AT" 

AHA ( I ) 
AB M I ) 
AR ? ( I ) 
AH 3 ( I } 


HK ( I ) *AR| ( I ) +GK ( I ) * A H ? ( I ) 

-GK ( I ) # A R 1 ( I ) 4 HK ( I ) * A : 3 2 ( l ) 

JK ( I ) * A R 1 ( I ) 4 < K ( I ) * A ■ J , 2 f ! ) + MK ( T ) * A 


UK ( I 5 * AB3 ( I ) -KK ( I ) * AB4 ( I ) 
K K ( [ ) * A B 3 ( I ) 4 JK ( I ) * A B 4 ( I ) 
3(1) +GK ( I ) * AB4 ( I ) 


KK ( I ) * AH 1 ( I ) - JK < I ) * A H? < 1 ) - 
AT 1 

A T 3 


( 1 ) * A H 3 < I ) + HK ( I ) * A 3 4 < ] ) 


C 

C* * x 


NORMAL I ZAT I ON 


1 


CORFCT 

I 

AH 1 ( I ) 
AB2 ( I ) 
AB 3(1) 
AB4 ( I ) 
RPT! )RN 
END 


! 


OF QUATERNION PARAMETERS 
O/SOPT ( A B 1 ( I ) * A R 1 ( I ) 4 Ar 


2 ( I } * A R ? ( I ) + AF-n ( I ) * AR3 ( I ) 


4 AR4 ( I ) * AR4 ( M ) 
CORECTt AR MI) 

C ORc. C T * AR2 ( I ) 

C 0 RtXT*AR 3 ( I ) 

C0RECT*AR4 ( I ) 
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APPENDIX B - Concluded 


c#*** 
c **** 

C**** 


SUBPOUT INF SAVE<I) 

SUBROUTINE FOR SAVING FORCING FUNCTION VALUES TO BE USED BY THE 
QUATERNION PARAMETER SUBROUTINE# REQUIRED WHEN A MULTI-PASS INTEGRATION 
ROUTINE IS USED IN ThE MAIN PROGRAM TO INSURE UPDATES ONLY ONCE 
PER ITERATION* 

COMMON/BK05/P fP)*Q(2)«R(2)* PDOT ( 2 > * QDOT ( 2 ) « ROOT ( 2 > 

COMMON/ TOU A T/PK ( 2 ) • Qk < 2 ) » RK < 2 ) ♦ PDOTK ( 2 ) • QDOTK ( 2 ) • RDOTK ( 2 ) 


PK ( ! >sP ( I ) 

SPDOTK ( 1 ) =POOT ( I ) 

OK ( I )=Q ( I ) 

SQDOTK ( I )=QOOT( T ) 

RK ( I ) = R ( I ) 

SRDOTKf I >=ROOT( I ) 

rftuon 


END 
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APPENDIX C 


A SIMPLIFIED ALGORITHM 

If one were to assume a zero-order hold on the p, q, and r inputs to the LL 
algorithm (that is, p, q, and r assumed constant over one iteration), a simplified 
version of LL may be obtained. This assumption implies that only the first two terms 
of equation (18) are used. As a result, the equations for A', B’, C', D 1 , Cj, and C 4 
would be eliminated. This simplification may be necessary for one with limited com- 
puter resources and would supply sufficient accuracy for low angular rates. The maxi- 
mum error on a design chart for this simplified algorithm (similiar to fig. 11) would be 
2.5 percent (or 9°). Typical errors are given in tables I, II, and III. It is important to 
notice that even this simplified LL algorithm is better than the classical second-order 
methods referred to in this report. 
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APPENDIX D 


STABILITY CHARACTERISTICS OF LL ALGORITHM 
The major computational difficulties associated with solving the system 

X = A(t) X(t) (Dl) 

where 

X(0) = x 0 

-r -q -p 

0 -p q 

p 0 r 

-q r 0 




X T X = 1.0 


using the classical integration methods in real time, lie in the fact that the eigenvalues of 
the system are outside the stability boundaries of the methods. The eigenvalues are 
purely imaginary with multiplicity two, that is, i, ± ^ i where = p2 + q2 + r 2 

Li U 

The basic problem is illustrated for AB-2 in figure 3 where the stability boundary (ref. 8) 
is shown relative to the eigenvalue location of equation (1). Clearly, the purely imaginary 
roots always lie outside the stability boundary of the method (h > 0) since the stability 
boundary only touches the imaginary axis at the origin. Therefore, any truncation or 
roundoff errors introduced in the solution process will, in general, not damp out; and the 
integral curve can become strongly divergent, increasing without bound as a function of 
time. All classical explicit methods which are suitable for real-time simulation have 
stability boundaries that do not include the imaginary axis except at the origin. 

The LL algorithm described previously can be written as 



(D3) 
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APPENDIX D - Continued 


where 


1 

0 

0 

0 


0 

G 

-J 

-K 

0 

1 

0 

0 

4- 

-G 

0 

-K 

J 

0 

0 

1 

0 


J 

K 

0 

G 

0 

0 

0 

1 


K 

-J 

-G 

o i 

i 


(D4) 


The method has local truncation error of 0(h 2 ). The scalars H, G, J, and K are as 
defined in equations (20) to (23). The physical meaning of X (ref. 3), indicates bounded 
solutions of the quaternion rate equations (neutrally stable solutions). Therefore, the 
eigenvalues of should lie close to the unit circle (ref. 9). Indeed a "perfect" algo- 

rithm would necessarily have root locations on the unit circle. The eigenvalues of the 
matrix can be written as 


?! = H + i y 

?2 = H - i y 

} 

? 3 = H + i y 

? 4 = H - iy 
^ J 


(D5) 


where 


7 - +\ 


G 2 + J 2 + K 2 


The magnitude of (j = 1, . . ., 4) is given as 




1) P(Z) + ®Q(Z) 


Jk 


where 


Y = 


to 


Z = 


k 11 


(D6) 


(D7) 

(D8) 
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APPENDIX D - Continued 


P(Z) = 


- + — (1 - cos Z) -MS_Z + lfez\ 
2 2Z 2 Z 4^ Z j 


(D9) 


Q(Z) = fe£ _ cos z 


(DIO) 


= ^:(PkPk + ^k + r k r k 


(Dll) 


The error in 


C« 


^ is defined as 


| ) P(Z) + (jW) 


(D12) 


Note for u> k * 0 and co k * 0, 


lim 

h~0 


f) p(z) + (f 1Q (Z ) 


= o 


J k 


(D13) 


that is, the method is consistent. Observe also for Y * 0, 


lim 

co— °° 



1.0 


(D14) 


when h is fixed and cl) § M (where M is a positive number), which is a desirable 
asymptotic property. For sufficiently small h 


cl) k o; k h J 1 . 

= 1 + + CO, 

6 16 k 


2 h 4 + o(h 5 ) 


(D15) 


with a corresponding error defined by equation (12). For the case where co is a constant, 
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APPENDIX D - Continued 


q 2 = 1.0 (D16) 

k 

indicating perfect placement of the system eigenvalues. Corresponding results for AB-2 
(ref. 1) and RK-2, two popular methods employed in real-time simulation programing, 
are, respectively, 



.4, ,4 
n /-v 

+ — - — + O 
32 



(D17) 




+ O 



(D18) 


For this case (c o is a constant) the superiority of the LL algorithm is obvious since there 
is a definite bias of root placement for the other two methods. The implication of this 
observation with regard to numerical stability will now be considered. 

Let X^tj^j be the exact solution vector for equation (1) at time t^ and be 

the one-step solution of the LL algorithm, and define the error as = X^ - X(t^j. It 
can be shown (ref. 10) that error propagation satisfies an equation of the form 


e k + l = M k e k + >fk’ h ) 


(D19) 


where y^tj c ,hj is the local truncation error. The behavior of for large n depends 
mainly on the eigenvalues of M^. If these eigenvalues have magnitudes which are con- 
sistently less than one, then the error norm will remain bounded. If u; is a constant 
and assuming no computer round-off errors, equation (D19) has the solution gj = 0 for 
all j. This conclusion is not true for AB-2 and RK-2 since the magnitude of the eigen- 
values are always greater than 1 for oo * 0. The analysis of equation (D19) in the general 
case is more difficult and will not be presented here. Experience has indicated, however, 
that the LL algorithm exhibits a high degree of stability for relative large step sizes, i.e., 
h = over a wide range of o> and cl). A possible explanation of this can be deduced 

U Li 

from equation (D15). In aerospace applications the acceleration must eventually 
change sign; therefore the eigenvalues of M k will not consistently lie outside the unit 

i 2 

circle. This error will accelerate when k . > 1 and will eventually damp when cl) 

changes sign. 
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APPENDIX D - Concluded 

Another indicator of performance of a specific numerical method with regard to the 
solution of equation (Dl) is the magnitude of the variation of the constant of motion X X. 
Defining V = X^X from equation (D3) and substituting equation (D4) there results in 

V k+1 = Cj 2 V k (D20) 

k 

It can be shown that V is a Liapunov function (ref. 11) for the system and therefore can 
be employed to study the stability characteristics of the LL algorithm, i.e., stability is 
implied by the condition (assuming no rounding errors) 

Q(Z) ^ 0 (D21) 

-k 

for all k. 
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1/32 sec unless otherwise noted. 




TABLE H.- TYPICAL ERRORS IN QUATERNIONS, DIRECTION COSINES, 
AND ATTITUDE ANGLES FOR SINUSOIDAL PULSE INPUTS 
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1/32 sec unless otherwise noted. 



TABLE III.- TYPICAL ERRORS IN QUATERNIONS, DIRECTION COSINES, 
AND ATTITUDE ANGLES FOR TAPED PILOTED RUN 
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h = 1/32 sec unless otherwise indicated. 




Euler angle to quaternion initialization 
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Figure 1.- Block diagram for determining aircraft attitude employing quaternion rate equations. 







Figure 2.- Block diagram comparing the implementation of AB-2 to that of LL for "norm in 
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(b) q = r = 0 and p = 1,3, and 5 rad/sec; h = 1/32 sec. 

Percent error in the norm squared of the Euler parameters using AB-2. 
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Time, t, sec 
(b) "Norm in." 

Figure 5.- Time-history segment of the quaternion aj, comparing AB-2 with LL 
for h = 1/32 sec and p = 10 rad/sec (q = r = 0). 
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Figure 8.- Differences in solution between AB-2 and LL for direction cosines, 
h = 1/32 sec; p = 10 rad/sec (q = r = 0); "norm in. M 
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(a) Time histories of the quaternions aj, . . a^ for AB-2 and LL. (The plus 
symbol denotes solution of RK-7 used as independent check.) 

Figure 9.- Sinusoidal inputs for h = 1/32 sec and "norm in." 
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(b) Time histories of the Euler angles <p, 6, and ip for AB-2 and LL. (The 

plus symbol denotes solution of RK-7 used as independent check.) 

Figure 9.- Continued. 
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(c) Differences in solution between AB-2 and LL for direction cosines 

Figure 9.- Continued. 
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(a) Time histories of the quaternions a^, . . a^ for AB-2 and LL. (The plus 
symbol denotes solution of RK-7 used as independent check.) 

Figure 10.- Sinusoidal pulse inputs for h = 1/32 sec and "norm in." 
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(b) Time histories of the Euler angles <f>, 9, and \p for AB-2 and LL. (The plus 

symbol denotes solution of RK-7 used as independent check.) 

Figure 10.- Continued. 
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Time, t, sec 

(c) Differences in solution between AB-2 and LL for direction cosines. 

Figure 10.- Continued. 
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(d) Differences in solution between AB-2 and LL for Euler angles. 

Figure 10.- Concluded. 
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(a) Direction cosines. 

Figure 13.- Differences in solution between AB-2 and LL for taped piloted 
run. h = 1/32 sec and "norm in.” 
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Figure 13.- Concluded. 
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(a) Time histories of the quaternions a 1? . . a 4 for AB-2 and LL. (The plus symbol 
denotes solution of RK-7 used as independent check.) 

Figure 14.- Sinusoidal inputs for h = 1/16 sec and "norm in.” 
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(b) Time histories of the Euler angles cp, 6, and if/ for AB-2 and LL. (The plus 
symbol denotes solution of RK-7 used as independent check.) 

Figure 14.- Concluded. 
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